rm(list=ls())
setwd("/Users/yuehou/Dropbox/KingReplication/Ongoing Work")
 

#############################
#############################
#### Loop with Beck-Katz ####
#############################
#############################

load("imputations23.Rdata")
m <- 8

# Loads initial imputations into dataframe
for (j in 1:m){
  assign(paste("data_",j,sep=""), amelia.out9$imputation[[j]])
}


for (j in 1:m){
	data_temp <- get((paste("data_",j,sep="")))
	

	# Reframe as a panel dataset
	library(plm)
	p_newset <- pdata.frame(data_temp)
	
	p_newset2<-subset(data_temp,dhasminwage==1)  #just for minwage
	p_newset2<-pdata.frame(p_newset2)
  #########Subsetting data only to include minwage=1
    
    p_newset2$log_oecd_gdp <- log(p_newset2$oecd_gdp)
    p_newset2$lag_ia5010 <- lag(p_newset2$final_ia5010)
	p_newset2$lag_govem_yh2<- lag(p_newset2$govem_yh2)

	p_newset2$lag_leftc<- lag(p_newset2$leftc)
	p_newset2$lag_leftc_2<- lag(p_newset2$leftc,2)
	p_newset2$lag_leftc_3<- lag(p_newset2$leftc,3)
	p_newset2$lag_leftc_4<- lag(p_newset2$leftc,4)
	p_newset2$lag_leftc_5<- lag(p_newset2$leftc,5)

	p_newset2$lag_leftgs   <- lag(p_newset2$leftgs)
	p_newset2$lag_leftgs_2 <- lag(p_newset2$leftgs,2)
	p_newset2$lag_leftgs_3 <- lag(p_newset2$leftgs,3)
	p_newset2$lag_leftgs_4 <- lag(p_newset2$leftgs,4)
	p_newset2$lag_leftgs_5 <- lag(p_newset2$leftgs,5)

	p_newset2$lag_arm_generos = lag(p_newset2$arm_generos)
	p_newset2$lag_minwag_jf = lag(p_newset2$minwag_jf)
	


# Transform generosity measure

	temp1 <- p_newset2$soc_exp1 <- (p_newset2$arm_generos - p_newset2$lag_arm_generos)
	temp2 <- p_newset2$lag_arm_generos
	temp3 <- p_newset2$soc_exp <- 100 * (temp1 / temp2)
	

	# Oil covariate
	p_newset2$oilshock <- 0
	
	p_newset2$oilshock <- (factor(p_newset2$year)==1973 | factor(p_newset2$year) == 1972 | factor(p_newset2$year) == 1971 | 	factor(p_newset2$year) == 1970)
	
 ###############################################
   

    # import model
   
   results <- plm(model="within",data=p_newset2,minwag_jf ~ lag_minwag_jf + vk_corp + vk_corp:leftgs + leftgs + vk_corp:lag_leftgs + lag_leftgs +  vk_corp:lag_leftgs_2 + lag_leftgs_2 + oecd_gdpgr + arm_unemp + arm_open + arm_finop + oecd_debt + oilshock )   #lag left 2
   
   
  results2 <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))
 
  # marginal effects (left and corporatism)
   hkvalue<-seq(1,5,1)

   beta_interaction<- results$coeff["leftgs"]+hkvalue*results$coeff["vk_corp:leftgs"]
   
   covariance1 <-vcovBK(results,cluster=c("group","time"))["leftgs","vk_corp:leftgs"] 
   
   
   se_interaction <- sqrt(results2["leftgs",2]^2+(hkvalue^2)*results2["vk_corp:leftgs",2]^2+2*hkvalue*covariance1)
   
	
	results <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))

	# Output vectors
	assign(paste("resultsc_",j,sep=""),results[,1])
	assign(paste("resultss_",j,sep=""),results[,2])
	assign(paste("resultsinteraction_",j,sep=""),beta_interaction)
	
	assign(paste("resultsinterSE_",j,sep=""),se_interaction)
}




# Vectors of Results   
coef <- 0
coef_interaction<-0

se <- 0
se2 <-0
se_interaction<-0
se_interaction_2<-0




for (j in 1:m){
	coef <- coef + get(paste("resultsc_",j,sep=""))
	coef_interaction<- coef_interaction + get(paste("resultsinteraction_",j,sep=""))
	se <- se + get(paste("resultss_",j,sep="")) 
	se2 <- se2 + (get(paste("resultss_",j,sep="")))^2
	se_interaction<- se_interaction + get(paste("resultsinterSE_",j,sep="")) 
	se_interaction_2<- se_interaction_2 + (get(paste("resultsinterSE_",j,sep="")))^2

}

coef <- coef/m
coef_interaction<-coef_interaction/m
se <- se/m
se2 <- se2/m
se_interaction_2 <-se_interaction_2/m
se3 <- 0
se4 <- 0


for (j in 1:m){
	se3 <- se3 + (get(paste("resultsc_",j,sep="")) - coef)^2
	se4 <- se4 + (get(paste("resultsinteraction_",j,sep=""))-coef_interaction)^2
}


se3 <- se3/(m-1)
se4 <- se4/(m-1)
final_se <- sqrt(se2 + se3*(1+(1/m)))
t <- coef/final_se

final_se_interaction<-sqrt(se_interaction_2 + se4*(1+(1/m)))

t_interaction<-coef_interaction/final_se_interaction

# Output Results
cbind(coef,final_se,t)


cbind(coef_interaction,final_se_interaction,t_interaction)



##### PLOT#####

vkrange=c(1,2,3,4,5)


plot(vkrange, coef_interaction,"l",lwd=2,ylim=c(-0.005,0.005),xlab="Corporatism",ylab="Marg. Effect of Left Gov",main="Minimum Wage")

lines(lowess(vkrange,coef_interaction + 1.96*final_se_interaction),col="red",lwd=1,lty=3)

lines(lowess(vkrange,coef_interaction - 1.96*final_se_interaction),col="red",lwd=1,lty=3)

abline(h=0,col="blue",lwd=.5)



